An analytic solution to LO coupled DGLAP evolution equations: a new pQCD tool 
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We have analytically solved the LO pQCD (leading order perturbative QCD) singlet DGLAP 
(Dokshitzer, Gribov, Lipatov, Alterelli, Parisi) equations P43l using Laplace transform techniques. 
Newly-developed highly accurate numerical inverse Laplace transform algorithms 0, H| allow us to 
write fully decoupled solutions for the singlet structure function F s (x,Q 2 ) and G(x,Q 2 ) as 

F 8 (x,Q 2 )=T s (F s0 (x ),G (x )) and G(x,Q 2 ) = G(F s0 (x ),G (x Q )), 

where the xo are the Bjorken-s: values at Qq. Here T s and Q are known functions — found us- 
ing LO DGLAP splitting functions — of the initial boundary conditions F s o(x) = F s (x,Qq) and 
Go(x) = G(x,Qo), i.e., the chosen starting functions at the virtuality Qq. For both G(x) and 
F 3 (x), we are able to either devolve or evolve each separately and rapidly, with very high numerical 
accuracy, a computational fractional precision of O(10 -9 ). Armed with this powerful new tool in 
the pQCD arsenal, we compare our numerical results from the above equations with the published 
MSTW2008 and CTEQ6L LO gluon and singlet F s distributions ffj, starting from their initial val- 
ues at Qq — 1 GeV 2 and 1.69 GeV 2 , respectively, using their choices of a s (Q 2 ). This allows an 
important independent check on the accuracies of their evolution codes and therefore the compu- 
tational accuracies of their published parton distributions. Our method completely decouples the 
two LO distributions, at the same time guaranteeing that both G and F s satisfy the singlet coupled 
DGLAP equations. It also allows one to easily obtain the effects of the starting functions on the 
evolved gluon and singlet structure functions, as functions of both Q 2 and Qq, being equally accu- 
rate in devolution (Q 2 < Q 2 ) as in evolution (Q 2 > Qo). Further, it can also be used for non-singlet 
distributions, thus giving LO analytic solutions for individual quark and gluon distributions at a 
given x and Q 2 , rather than the numerical solutions of the coupled integral-differential equations 
on a large, but fixed, two-dimensional grid that are currently available. 



I. INTRODUCTION 



The search for new physics at the LHC demands an accurate knowledge of gluon distribution functions at small 
Bjorken x and large virtuality Q 2 , both for estimating QCD backgrounds and for calculating gluon-initiatcd processes. 
The traditional method has simultaneously determined gluon and quark distribution functions by fitting experimental 
data on neutral- and charged-current deep inelastic scattering processes and some jet data over a large domain of 
values of x and Q 2 . The distributions at small x and large Q 2 are determined mainly by the proton structure 
function K^x, Q 2 ) measured in deep inelastic ep (or 7*p) scattering. The fitting process starts with an initial Qq, 
typically less than or equal to the square of the c quark mass, m 2 « 2 GeV 2 , and individual quark and gluon trial 
distributions parameterized with pre-determined shapes, given as functions of x for the chosen Qq. The distributions 
are then evolved numerically on a finite, albeit large, two-dimensional grid in x and Q 2 to larger Q 2 using the coupled 
integral-differential DGLAP equations typically in leading order (LO) and next-to- leading order (NLO), and 

the results used to predict measured quantities. The final distributions are then determined by adjusting the input 
parameters to obtain a best fit to experimental data, fitting both HERA and Tcvatron data over a large range of 
x and Q 2 , along with selected hard scattering data from fixed target experiments. This procedure is very indirect 
in the case of the gluon: the gluon distribution G(x,Q 2 ) = xg(x,Q 2 ) does not contribute directly to the accurately 
determined structure function F^ p (x, Q 2 ), and is determined only through the quark distributions in conjunction with 
the evolution equations, or at large x, from jet data. For recent determinations of the gluon and quark distributions, 
see [l-Eil. 

In the following, we will summarize our method for analytically determining G{x, Q 2 ) and the singlet structure 
function F s (x,Q 2 ) directly and individually, using as input Gq(x) = G[x,Qq) and F s q{x) = F s (x,Qq), where Qq is 
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arbitrary, with the guarantee that each distribution individually satisfies the coupled DGLAP equations. The method 
is readily extended to embrace non-singlet functions, so that it can be used also to find individual quark distributions. 
However, we will not pursue that goal in this communication. Instead, we give a numerical demonstration which takes 
advantage of the fact that our analytic solutions achieve numerical accuracies of O(10 -9 ), giving us a new diagnostic 
tool to verify published LO singlet structure functions (F s (x, Q 2 )) and gluon (G(x, Q 2 ) = xg{x, Q 2 )) distributions. In 
order to test the numerical accuracy of their evolution codes, we consider two cases, using the published LO starting 
distributions for Go and F s q: 

1. MSTW2008 @: for Qq = 1 GeV 2 , we generate LO singlet structure functions and gluon distributions [f|, 
using the strong coupling constant a s (Q 2 ) that they used for their LO evolution and compare them with their 
published values MSTW2008 [<J for the domain 1CT 6 < x < 1 and 1 < Q 2 < 100000. We find that their 
evolution code has serious problems at small x, producing significant numerical inaccuracies. It should be noted 
that the MSTW group does not do devolution. 

2. CTEQ6L 0: for Q 2 = 1.69 GeV 2 j_we generate LO singlet structure functions and gluon distributions, using the 
strong coupling constant a s (Q 2 ) [8| they used for both LO evolution and devolution. With our high numerical 
precision at all x and Q 2 , we are able to verify all of their published evolution results — to larger Q 2 — but show 
that their published devolution results, i.e., Q 2 < Qq, have significant numerical inaccuracies at small x. 

Finally, using our accurate CTEQ devolution results, we compare LO starting distributions for both groups at Q 2 = 1 
GeV 2 , noting that the CTEQ6L LO gluon distribution turns over and goes negative at small x, i.e., x < 5 x 10 -5 , 
whereas the MSTW2008 LO gluon starting distribution continues to rise sharply at small x. 



II. DECOUPLING THE COUPLED LO SINGLET DGLAP EQUATIONS 

Our approach uses a somewhat unusual application of Laplace transforms (l2l . Il3| , in which we first introduce the 
variable v = ln(l/ic) into the coupled DGLAP equations, then Laplace transform these coupled integral-differential 
equations in v space to obtain coupled homogeneous first-order differential equations in the Laplace-space variable s. 
We solve these equations analytically. Finally, using fast and accurate numerical inverse Laplace transform algorithms 
we transform the solutions back into v space, and, finally, into Bjorken x-space, so that we can write 

F s {x,Q 2 )=T s (F s0 (x Q ),Go(x )) and G(x,Q 2 ) = Q(F s q(x ),Go(x )), (1) 



where the functions T and Q arc determined by the splitting functions in the DGLAP equations, with xq being the 
Bjorkcn-x at the starting virtuality Qq; F s q(x) and Gq(x) are the known starting distributions at Q 2 = Qq, where 
evolution (devolution) begins. 

Our method can be generalized to NLO (see Ref. [Hj]), but for brevity, we will limit ourselves to LO in this paper. 
We write the coupled LO DGLAP equations [HI, [Jj| as 

^ 9F s _ AV , ^ ,16 ^,.1-1, 16. f 1 (F s (z,Q 2 ) F s (x,Q 2 )\ dz 



-(x,Q 2 ) = AF s (x,Q 2 ) + -F s (x,Q 2 )\n + —x 

a s (Q z ) omQ z 6 x i J x \ z x 

-^^ 1 F > (3,<3 2 )(l + 2)^ + 2n / x £ G (;,0 2 ) (1-2^ + 2^) J, (2) 
a s (Q z ) olnQ z 3 x Jx \ z x ) z — x 

+-f ^> (; - 2 + ; - S) d i + If ^ (^(i- ) 2 ) t • < 3 » 

Here a s (Q 2 ) is the running strong coupling constant, and for LO MSTW2008 @ is given by the LO form 

with rif the number of quark flavors. The QCD parameter A5 is fixed so that the known a s (A/§) is reproduced and 
then A4 and A 3 are adjusted so that a s is continuous across the boundaries Q 2 = M 2 and M 2 , respectively, where 
Mb and M c are the masses of the b and c quarks. Later, we will also introduce the NLO form of a s used (with 
a s (M§) = 0.118) in their LO CTEQ6L evolution, when we discuss CTEQ6L pdfs. 
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We now examine the last two terms of line 1 in Eq. and rewrite them, introducing the variable changes v 
ln(l/x), w = ln(l/z), and the notation F s (v,Q 2 ) = F s (e~ v ,Q 2 ), G(v,Q 2 ) = G(e~ v ,Q 2 ), as 



-F s (v,Q 2 )\n(e v - 1) + y J {F s (w,Q 2 ) 
= y^|^KQ 2 )ln(l-e-(— )) dw. 



F s (v,Q A )e 



2 \ _tf— w 



dw 



(5) 



where the final result — the last line in Eq. ([5]) — is found by replacing the upper limit v in integral of line 1 of Eq. ([5]) 
by v — e, carrying out the integrals, doing a partial integration and finally, taking the limit as e — > 0. Similarly, we 
find for the last two terms of line 1 in Eq. that 



12G(v, Q 2 ) ln(e° - 1) + 12 / G(w, Q 2 ) - G(v, Q 2 )e 



\2 \ v—w 



dw 



= 12 



v dG 
o dw 



(w,g 2 )ln(l-e- (, '- u ' ) ) dw. 



(6) 



We now rewrite Eq. ([2]) and Eq. (j3]) in terms of the new variable v as 



4tt dF s 
a s (Q 2 )d\nQ 2 



(v,Q 2 ) = 4F s {v,Q 2 ) + 



16 r OF* 
dw 



{w,Q 2 )\n(l-e w - v ) dw 
(w,Q 2 ) (e-( v ~ w) + e" 2 ^-™)) dw 
+2n f x I G{w, Q 2 ) (e- (v ~ w) - 2e- 2(v ~ w) + 2er z{v - w) ^j dw, 



F„ 



(7) 



4tt dG 
a s (Q 2 )d\nQ- 



iv,Q 2 ) 



33 - 2n f 



dG 



G{v,Q 2 ) + 12 y _( w ,Q 2 )ln(^l-e- (l '-^j 
+ 12 f G(w, Q 2 ) (l - 2e- {v - w) + e - 2{v ~ w) - e - 3( - v ~ w ^ dw 



F s (w,Q 2 ) 1+ 1-e 



G?U>. 



(8) 



The DGLAP equations have now been written in a form such that all of the integrals in Eq. (J7]) and Eq. (|SJ) arc 
manifestly seen to be convolution integrals. Thus, introducing Laplace transforms allows us to factor these convolution 
integrals, since the Laplace transform of a convolution is the product of the Laplace transforms of the factors, i.e., 



w] dw, s 


= C 


[L 









r 

/ F[v — w]H[w] dw; s 
Jo 



C 



Defining the Laplace transforms of F s (v,Q 2 ) and G(v,Q 2 ) in s space as 



C[F[v];s]xC[H[v];s]. 



(9) 



f(s,Q 2 ) ee C[F s (v,Q 2 ); S 
and noting that 



F s (v, Q 2 )e~ sv dv, g{s, Q z ) = C[G{v, Q 2 );s] 



dw 



(w,Q 2 );s 



sf(s,Q 2 ), 



C 



dG 

dw 



(w,Q 2 );s 



= sg{s,Q 2 ), 



G(v,Q 2 )e~ sv dv (10) 



(11) 



since F s (v = 0, Q 2 ) = G(v = 0, Q 2 ) = 0, we now factor the Laplace transforms of Eq. (0 and Eq. (jHJ into two coupled 
first order differential equations in Laplace space s having Q 2 -depcndent coefficients. These can be written as 



df 



dlnQ- 



is,Q 2 ) 



a s (Q 2 

47T 



-$ / ( s )/( s ,Q 2 ) + 



a s (Q 2 ) 



-in 



@ f (s)g(s,Q 2 ), 



dg 

d\nQ' 



is,Q 2 ) 



4tt 



-* g (s)g(s, Q 2 ) + -J^Q g ( s )f(s, Q 2 ). 



(12) 
(13) 
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The coefficient functions $ and 9 arc given by 



*/(*) 

e f ( s ) 



4 - 
2» 



1 



1 



8 

3 \s+l ' s+ 
1 2 



33 - 2n 



+ 2 + 1) + 7b) 
2 



1 s + 2 s + 3 



/ 



3 

8 /2 



12 



1 



1 



3 V s 



1 2 
s s + 1 + s + 2 s + 3 
1 



V>(s + 1) - 7b 



1 s + 2 



(14) 
(15) 
(16) 
(17) 



where ip(x) is the digamma function and = 0.5772156 ... is Euler's constant. 

The solution of the coupled equations in Eq. (fT2j) and Eq. (|T3| in terms of initial values of the functions / and 
<7, specified as functions of s at virtuality Qq, is straightforward. The Q 2 dependence of the solutions is expressed 
entirely through the function 



r{Q\Ql) 



4tt 



Q 



a s (Q' 2 )d In Q' 



With the initial conditions fo{s) = f(s, Qq) and go{s) = <?(s, Qo)> the solutions are 

f(s, T ) = kff{s,T)fo(s) + kfg(s,T)g (s), 
g(s,r) = k gg (s,T)g (s) + fc s /(s, r)/ (s), 



where the coefficient functions in the solution are 

i(4/W+« s W) 



fc //( S , T ) 

kfg(s,T) 
kgg(s, T) 

k gf {s,T) 



cosh Q^oo) 

ei($/(s)+4 , s(s)) 2sinh(§ J R( S )) 



sinh (fi?(s)) 



R(s) 



i?( s ) 



©/(*) 



cosh (*,(-) -*.(.)) 



i?(.s) 



e 2 



2 sinh (§i?(s)) 



e s (s), 



(18) 



(19) 
(20) 



(21) 
(22) 
(23) 
(24) 



with i?(s) = J ($/(s) — $ 9 (s)) 2 + 40 f(s)@ g (s). Clearly, the fundamental solutions in Laplace space s, Eq. (fl~9|) and 
Eq. (|20[) . are symmetric under the interchange f <-> g. 

Let us now define four kernels Kpp, Kfg, K G p and K GG , the inverse Laplace transforms of the fc's, i.e., 

K ff (v,t) = L-^kff^T^v], K FG (v 1 T) = £- 1 [k f g(s 1 T);v] 1 (25) 
K gg (v,t) = L-^kgg&T^v}, *r OF (^r) t);v]. (26) 

It is evident from Eqs. (fl"8|) . (f22|). and ([24]) that Kp G and -FTgf vanish for Q 2 = Qq where t(Q 2 ,Qq) = 0. It can also 
be shown without difficulty that for r = 0, Kpp(v,0) = Kgg{v,0) = S(v) and that Kp G (v,0) = K G p(v,0) = 0. 

The initial boundary conditions at Qq are given by F s q(x) — F s (x,Qq) and Gq(x) = G(x, Qq). In w-space, 
F s0 (v) = F s0 (e~ v ) and Gq(v) = Go(e~ v ) are the inverse Laplace transforms of /o(s) and go(s), respectively, i.e., 



F s0 (v) = C-^Ms^v] andGo(w) = [ffo(s); 



(27) 



Finally, we can write our decoupled singlet structure function F s and G solutions in w-space in terms of the convo- 
lution integrals as 



Fs(v,Q 2 ) = 
G(v,Q 2 ) = 



(v-w,T(Q 2 ,Q 2 Q ))F s0 (w)dw+ / K FG (v-w,T(Q 2 ,Q 2 ))G {w)dw, 



L GG 



(v - w, t(Q 2 ,QI))G {w) dw + 



*-GF 



(v-w,T(Q 2 ,Q 2 ))F s0 (w)dw. 



(28) 
(29) 
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We now derive an alternate form of the solution to the decoupled equation, very useful for computational purposes, 
that does not use the convolution theorem. Using a suitable fast and accurate numerical inverse Laplace transform 
[H, we can directly invert Eq. (|T9|) and Eq. (|20|) . since /o(s) and go(s) — the Laplace transforms of the known starting 
functions F s0 (v,Q 2 ) and Gq(v,Q 2 ) — are readily obtainable; the coefficient functions, the fc's given in Eq . (f2"TH2"4")) . are 
known functions of s and r, and hence, of Q 2 and Qq. Thus we finally write our decoupled analytic solution in v 
space as 

F s (v,Q 2 ) = C- 1 [(kff(s,T)Ms) + k fg (s,T)go(s));v), (30) 
G{v,Q 2 ) = C- 1 [(k gg ( S ,T)go(s) + k gf ( S ,T)f (s));v}. (31) 

In order to use our solution in the integral representation of Eq. (|28|) and Eq. (|29[) , we must first numerically invert 
Laplace transforms of the type of kff and k gg that for small r look similar to Dirac 8 functions; a formidable numerical 
task that is inherently inaccurate, and is thus computationally intensive and significantly slower (but possible) using 
the numerical inverse transforms of Ref. [j| ■ On the other hand, if we use Eq. (|3T))) and Eq. pip , we only have to invert 
a function whose inverse Laplace transform (F s (v, Q 2 ) or G(v, Q 2 )) is very smooth and thus can be well approximated 
by a high order polynomial in v. As shown in Ref. |5|, it can then in principle be evaluated to arbitrary accuracy very 
rapidly. It will be shown in the Appendix that we actually achieve a fractional accuracy of O(10~ n ) in our numerical 
Laplace inversion. In Section EH we will do a detailed evaluation of the inherent overall numerical accuracy for actual 
physical problems, showing that we can do both devolution and evolution rapidly to fractional accuracies of O(10 -9 ) 
using the numerical methods outlined in the Appendix. 

The final desired decoupled F s (x, Q 2 ) and G(x, Q 2 ) in Bjorkcn-x space are readily found by substituting v — ln(l/a;) 
into the w-space solutions for F B (v, Q 2 ) and G(v, Q 2 ) from Eq. (|30|) and Eq. (f3T1) . 



III. ANALYTIC LO NON-SINGLET DISTRIBUTIONS 



For non-singlet distributions F ns (x,Q 2 ), such as the difference between the u and d quark distributions, 
x [u(x, Q 2 ) — d{x 1 Q 2 )] , we can schematically write the logarithmic derivative of F ns as the convolution of F ns (x, Q 2 ) 
with the non-singlet splitting function K. ns {x) (using the convolution symbol (§)), i.e., 



4tt 



dF n 



-(x,Q 2 ) 



F„ 



(32) 



a s (Q 2 )dHQ 2 ) 

After again changing to the variable v = ln(l/x) and going to Laplace space s, we find the simple solution 

/„ s (s,r) = e r *" s(s) / ns o(s), where $ ns (s) = C e~ v IC ns (v);s and K. ns (v) = K, ns (e~ v ) . (33) 

Thus we can find any non-singlet solution in w-space, using the non-singlet kernel K ns (v) = Cr 1 [e T * rea ( s ) ; v\ , by cither 
employing the Laplace convolution relation 



F ns {v,Q 2 ) 



K ns (v - w,T(Q 2 1 Ql))F ns0 (w)dw 



or the non-integral form 



F ns (v,Q 2 ) = £- 1 \e T<s, ^ s \f ns o( S );v 



(34) 



(35) 



In this case, either method works equally well numerically, since the non-singlet functions K ns (v) can also be approx- 
imated by a polynomial in v. 

For brevity, we will not pursue the case of the non-singlet solution any further here except to note that in LO the 
$„ s (s) in Eq. (|3"3")l is identical to &f(s) defined in Eq. (fT4")) . Instead, we will concentrate on the more difficult case of 
F s and G. 



IV. LO MSTW2008 SINGLET AND GLUON DISTRIBUTIONS 



As an example of the application of our analytic decoupled solutions, we will use the published MSTW2008 initial 
starting functions F s q(x) and Gq(x) at Qq = 1 GeV 2 Q and will compare our LO x-space gluon distribution G{x, Q 2 ) = 
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xg(x, Q 2 ) using Eq. (|31l) and our LO singlet structure function F s (x, Q 2 ) using Eq. (|30|) — both numerically evaluated 
using a powerful new inverse Laplace transformation algorithm Q — with the corresponding LO distributions published 
by the MSTW collaboration Q. In order to insure continuity across the boundaries Q 2 = M 2 and M 2 , we will first 
evolve from Ql = 1 GcV 2 (the MSTW Q 2 value) to M 2 and use our evolved values of G(x, M 2 ) and F s (x,M 2 ) for 
new starting values Gq(x) and F s q(x). We will then evolve to M 2 , repeating the process, thus insuring continuity 
of G and F s at the boundaries where nj changes. We use the MSTW values M c = 1.40 GeV, Mb = 4.75 GeV, 
a s (l GcV 2 ) = 0.6818 and a s {M 2 ) = 0.13939 in their definition of a s (Q 2 ) in Eq. (®. 



A. G(x,Q 2 ) and Fs{x,Q 2 ) for LO MSTW2008 

In Fig. [T]we show the LO x-space results for G(x,Q 2 ) = xg(x,Q 2 ) (upper figure) and F s (x,Q 2 ) (lower figure) vs. 
x, for 4 representative values of Q 2 . The x-domain, 10 -6 < x < 1, is the complete region covered by the MSTW 
group Q. The curves are the published LO MSTW2008 distributions @: from bottom to top; the (red) curve is for 
Q 2 = 5 GeV 2 ; the (brown) dashed curve is for Q 2 = 20 GeV 2 ; the (blue) dot-dashed curve is for Q 2 = 100 GeV 2 ; 
the (black) dotted curve is for Q 2 = M§. The (red) dots are our analytic results for LO G(x,Q 2 ) from Eq. (fBTj) and 
F s (x,Q 2 ) from Eq. (|3"0"j) . converted to x-space, using the LO MSTW2008 values for F s0 (x) and Go(x); the numerical 
values were evaluated using Mathematica [l5j . An outline of the numerical procedure is given in the Appendix. 

For large x, the agreement is excellent for all Q 2 . However, as seen in a close inspection of Fig. [I] the disagreement 
for both G and F s becomes significantly large as we go to small x. We will explore this in detail in Section TlV Bl 
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FIG. 1: Plots for LO MSTW2008 gluon distributions G(x,Q 2 ) = xg(x,Q 2 ) (upper plot) and F s {x,Q 2 ) distributions (lower 
plot) vs. Bjorken x. The MSTW2008 curves are for Q 2 = 5, 20, 100 and M% GeV 2 , bottom to top. The (red) dots are our 
evolution results for LO G(x,Q 2 ) from Eq. (J3TJ and F s from Eq. ([50)) . after converting to x-space, using the LO MSTW2008 
[f| values for F a o(x) and Go(x), with their choice of Q\ — 1 GeV 2 . The x range covers all of the published LO MSTW2008 x 
data. 
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B. Accuracy of evolved LO MSTW2008 distributions 



We now investigate quantitatively the accuracy of the evolved LO MSTW2008 distributions (Q 2 > Qq), introducing 
the fractional accuracy variable 

Fractional Accuracy = 1 - /^bdhm/Amstw, i = 1,2, (36) 

where fi = F s , f% = G. with BDHM denoting our LO analytic evaluations and MSTW denoting the published 
LO MSTW2008 values [f|. We show in Fig. [2] the fractional accuracy for the LO MSTW published distributions 
@ G(x,Q 2 ) (upper figure) and F s (x,Q 2 ) (lower figure) using the same four Q 2 values and legends used in Section 
HVl and Fig. [TJ i.e., the (red) curves are Q 2 = 5 GeV 2 ; the (brown) dashed curves are Q 2 = 20 GeV 2 ; the (blue) 
dot-dashed curves are Q 2 = 100 GeV 2 ; the (black) dotted curves are M§. Both the MSTW2008 G and F s are in 
excellent agreement with our (much more numerically precise) calculations in the domain x > 10 , with a fractional 
accuracy of ~ 0.1 — 0.5%. However, as is clearly seen in Fig. [T]for both G and F s and for all Q 2 , there is the same 
inaccuracy pattern in x, an increase of the fractional accuracy to ~ 2% down to x « 8 x 10 -6 , followed by a dip at 
x 4 x 10 -6 , with a hnal rise to another maximum at x w 2 x 10 -6 whose fractional accuracy is ~ 12%. These 
final inaccuracies at small x are quite significant. Since the x patterns are essentially independent of whether we 
are evaluating either G or F s , as well as being independent of Q 2 , they suggest that the MSTW numerical program 
undergoes a significant structural change at some unique value of x, independent of Q 2 , that seriously degrades their 
numerical output, leading to large errors at small x. The largest errors occur at the smallest Q 2 ; at Q 2 , (not shown) 
the error is ~ 12 — 13 %, and decreases monotonically to ~ 4 — 5 % at the highest Q 2 . As we will later see in Section 
[VT there is no such pattern in the LO CTEQ6L data 0j- 




FIG. 2: Fractional accuracy plots for LO MSTW [f| gluon distributions G(x,Q 2 ) = xg(x,Q 2 ) (upper plot) and F a (x,Q 2 ) 
distributions (lower plot), for Q 2 — 5, 20, 100 and M| GeV 2 , where the fractional accuracy is given by Eq. I|36p . The (red) 
curves are Q 2 = 5 GeV 2 ; the (brown) dashed curves are Q 2 = 20 GeV 2 ; the (blue) dot-dashed curves are Q 2 = 100 GeV 2 ; the 
(black) dotted curves are M 2 Z . The x range covers all of the published LO MSTW x data. 
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V. LO CTEQ6L SINGLET AND GLUON DISTRIBUTIONS 

As a second example of the application of our analytic decoupled solutions, we will compare our LO x-space gluon 
distribution G(x,Q 2 ) = xg(x,Q 2 ) from Eq. (|31|) and our LO singlet distribution function F s (x,Q 2 ) from Eq. (|3TJf — 
using the published LO CTEQ6L Q initial conditions at Qq = 1.69 GeV 2 — with the corresponding LO CTEQ6L 
distributions Q. In order to insure continuity across the boundary M 2 , we will first evolve from Qq = 1.69 GeV 2 
(the CTEQ6L Qq value) to M 2 and use our evolved values of G(x, M 2 ) and F s (x, M 2 ) for new starting values Go (a;) 
and F s q(x), thus insuring continuity of G and F s at the boundary where rif changes. We use the CTEQ6L values 
M c = 1.3 GeV and Mb = 4.5 GeV. We here use a NLO version of a s (Q 2 ), with a s (Af§) = 0.118, made continuous at 
Mb and M c , that was utilized in CTEQ6L (for details see Ref. Q). 




10~ 6 10~ 5 l(T 4 0.001 0.01 0.1 



x 



FIG. 3: Plots for LO CTEQ6L gluon distributions G(x,Q 2 ) = xg{x,Q 2 ) (upper plot) and F s (x,Q 2 ) distributions (lower 
plot) vs. Bjorken x. The curves are for Q 2 = 10 ,22, 90, 1200 and M 2 Z GeV 2 , bottom to top. The (red) dot s are our evolution 
results for LO G(x,Q 2 ) from Eq. f|31|) and F s from Eq. (|30[) (converted to x-space) using the LO CTEQ6L values for F B o(x) 
and G (x), where Q% = 1.69 GeV 2 . The x range in this Figure covers all of the published LO CTEQ6L x data. 



A. G(x,Q 2 ) and Fs(x,Q 2 ) for LO CTEQ6L 

In Fig. [3]we show the Bjorken x-space results for LO G(x, Q 2 ) = xg(x, Q 2 ) (upper figure) and LO F s {x, Q 2 ) (lower 
figure) vs. x, for 5 representative values of Q 2 . The x-domain, 1CV 6 < x < 1, is the complete region covered by the 
CTEQ group Q. The curves are the published CTEQ6L Q LO distributions. From bottom to top; the (red) curve 
is for Q 2 = 10 GeV 2 ; the (brown) dashed curve is for Q 2 = 22 GeV 2 ; the (blue) dot-dashed curve is for Q 2 = 90 
GeV 2 ; the (black) dotted curve is for Q 2 = 1200 GeV 2 ; the (orange) curve is for Q 2 = M§. Since CTEQ6L started 
evolution at Qq = 1.69 GeV 2 , we used F s0 and Go constructed from their values at Qq = 1.69 GeV 2 in Eq. (pIT]) 



9 



and Eq. ([30]). The (red) dots are our results for LO G(x, Q 2 ) from Eq. ([31]) and F s (x, Q 2 ) from Eq. ([30]) converted to 
x-space, using LO CTEQ6L values for F s0 (x) and Go(x), evaluated using Mathematica |ZL5j | . 

For all Q 2 the agreement is excellent over the entire x region, with a fractional accuracy of about ±5 x 10 -4 , 
(completely consistent with the 4 significant figures that are published) — for all F s and G at the five virtualities that 
we evaluated — with a minor and numerically unimportant exception of the lowest x region of F 8 (x, Q 2 = 22), where 
there was an offset of w 2 x 10~ 3 . 

B. Accuracy of CTEQ6L devolved distributions 

In Fig. [3j all of the distributions were for evolutions of G and F s from the CTEQ6L Ql = 1.69 GeV 2 to larger 
Q 2 . For another physics investigation, not relevant to this paper, we decided to compare LO starting distributions for 
MRSTW2008 and CTEQ6L at the MSTW2008 starting value of Ql = 1 GeV 2 . Using n f = 3, we devolved G and F s 
from the CTEQ6L starting values at Ql = 1.69 GeV 2 down to Q 2 = 1 GeV 2 , the MSTW2008 starting value for Ql. 

The results of this devolution are compared to the published CTEQ6L values Q in Fig. @] for G (upper figure) 
and F s (lower figure). In all cases, when we refer to "published CTEQ6L values", we mean the results found on 
the Durham pdf generator web site; see footnote [I?} ■ The solid (black) curves arc for CTEQ6L and the (red) dots 
are from Eq. (pTTj) and Eq. ([3"U| . In marked contrast to their evolution results, the CTEQ6L devolution results are 
numerically unstable, with F s being wrong by « 12% at x = 10~ 6 . We also note that there are large disagreements 
with their devolved G(x) for small x. Clearly, they have chosen to chop off their G distribution at small x, i.e., to write 
G{x) = for small x, rather than allow it to become negative. The errors for both G and F s become insignificant 
as x approaches 1. It is clear that CTEQ encounters major problems with the numerical stability of their published 
results for Q 2 < Qq, whereas they are completely accurate for Q 2 > Ql- 

For comparison, we also show in Fig. @] the published MSTW2008 starting distributions @ Go(x) and F s q(x) at 
Ql = 1 GeV 2 , the dashed (blue) curves. We note that the LO gluon distributions of the two different collaborations, 
when evaluated at the same virtuality, Q 2 = 1 GeV 2 , bear little or no resemblance to each other, with the CTEQ6L 
gluon distribution going negative for x < 3 X 10~ 5 . Although both singlet structure functions F s (x 7 Q 2 = 1) stay 
positive — as they must — Fig. 0] shows that there are also large differences between the two singlet structure functions 
at low x. 



VI. OVERALL NUMERICAL ACCURACY OF ANALYTICAL DEVOLUTION AND EVOLUTION 

As mentioned in Section IV B[ we had devolved from Ql = 1.69 GeV 2 to Q 2 = 1 GeV 2 , using the known CTEQ6L 
Go(x) and F s q(x) starting values. To estimate the overall accuracy of our entire numerical procedure, we took our 
devolved distributions G(x, Q 2 = 1) and F s {x 1 Q 2 = 1) and used them as starting values so that we could again evolve 
back to Q 2 = 1.69 GeV 2 . Finally, we compared the evolved numerical results with the original F s0 (x) and Gq{x), the 
distributions that we started with at Q 2 = 1.69 GeV 2 . An outline of our entire numerical procedure is given in the 
Appendix. 

In Fig. [5J we show the fractional accuracy of this "round-trip" comparison. The upper figure is for G and the lower 
figure is for F s . The (red) dots are the "round-trip" fractional accuracies at discrete x- values chosen to start and end 
this numerical exercise (corresponding to the transformed zeroes of the Chebyshev polynomials that we discuss in the 
Appendix). For the visual convenience of the reader, we have connected the dots. 

Where either G and F s is significantly large (x < 0.3), we see that the "round-trip" error is <4x 1CU 9 , thus 
yielding an overall error estimate of < ± 2 x 1CU 9 for either evolution or devolution. Detailed causes for this error 
are discussed in the Appendix. 

It is gratifying that the overall numerical uncertainty in our LO analytically decoupled solutions is small, thus 
furnishing us not only with a new accurate and fast calculation tool for exploring the effects of the shapes of different 
starting value distributions, but also with a diagnostic tool for easily determining the numerical calculational reliability 
of the already published parton distribution functions that are currently in major use by the high energy physics 
community. 

VII. CONCLUSIONS 

In conclusion, we have constructed decoupled analytical solutions for F s (x,Q 2 ) and G(x,Q 2 ) from the coupled LO 
DGLAP equations, yielding accurate numerical results for both evolution and devolution of O(10 -9 ) — a fast tool to 
study the dependence on the shape of the starting distributions -F^o^) and Gq(x), the boundary conditions at the 
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FIG. 4: Plots of gluon distributions G(x,Q 2 ) = xg(x,Q 2 ) (upper plot) and F s (x,Q 2 ) distributions (lower plot) vs. Bjorken 
x, at the devolved value of Q 2 = 1 GeV 2 . The (red) dots are our devolution results; the (black) solid curves are the published 
CTEQ6L results @] and the (blue) dashed curves are the starting F s0 {x) and Go (a;) for MSTW2008 @]. 



starting value Qq. Similar procedures can be used for non-singlet distributions, allowing one to obtain analytic LO 
solutions for individual quark distributions, as well as for the gluon distribution; thus avoiding the necessity for purely 
numerical solutions of the coupled DGLAP equations on a giant two-dimensional grid in (x, Q 2 ) space. In essence, 
using a program such as Mathematica [l5j], we can now define a parton distribution function for each quark and 
gluon and — after inputting the desired x and Q 2 — evaluate it accurately and rapidly (for a fast Mathematica program 
calculating LO F s (x,Q 2 ) and G(x,Q 2 ), see the Appendix). 

We have also used our analytic solutions coupled with the MSTW2008 initial starting functions [f| as a new and 
powerful diagnostic tool to study the numerical accuracy (the computational accuracy of their evolution code) of the 
LO MSTW2008 published distributions 0. For the small x-region, x < 10~ 4 , we discovered a pattern of significant 
numerical (computational) errors for both F s and G, ranging up to ~ 12% at the smallest x values in the published 
MSTW2008 results ®, true for all Q 2 . 

Applying the same new tools to CTEQ6L, we found no errors (to their accuracy of 4 significant figures) in either 
F s or G values when they did evolution from Qq = 1.69 GeV 2 to higher Q 2 values, but significant errors — increasing 
with decreasing x — when they did devolution to smaller Q 2 . In the future, we intend to evaluate F s q(x) and Gq(x) 
in both LO and NLO, from a fit to small x experimental data for the structure function F^^x, Q 2 ), in order to 
obtain (analytically) accurate values of G(x, Q 2 ) directly tied to experiment, which are needed for the interpretation 
of experiments at the LHC. 
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FIG. 5: Fractional accuracy plots for our LO gluon distributions G(x,Q 2 ) — xg(x,Q 2 ) (upper figure) from Eq. (|31[) and 
F s (x,Q 2 ) distributions (lower figure) from Eq. (|30|) . These accuracy estimates resulted from devolution from Q 2 = 1.69 GeV 2 



to Q A = 1 GeV 2 , then using these results for evolution back to Q 2 
from comparing the original values with the devolved-evolved ones. 



GeV . The fractional value error estimates result 
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Appendix A 

We outline here the actual calculation procedures necessary for fast and accurate numerical evaluations of Eq. ([50)1 
and Eq. pip . These calculations, although robust, require delicate choices as to the numerical techniques used in 
evaluating Eq. (J30J) and Eq. pTj) . 

As shown in Ref. if the function g(s) goes to at oo more rapidly than 1/s, then we can accurately approximate 
its inverse Laplace transform G(v) by 

2 N 

G(v) ~ -VRe^^)], (Al) 
v z — ' 

i=l 

where 2N is the order of the approximation, Wj and a,;, i = 1,2, ... , 2N, are known complex numbers for a given 2N, 
occurring in complex conjugate pairs. The actual numerical evaluation of Eq. (|A1[) can be quite unstable if one doesn't 
utilize arbitrary accuracy arithmetic as discussed in Ref. Q, since the weight functions u>i become exceedingly large, 
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even for modest 2N, and oscillate in sign Q. The use of Mathematica (or similar programs, which also carry out 
arithmetical operations to arbitrary accuracy) makes this requirement easy to satisfy. 

As shown in Ref. the inverse Laplace transform approximation to G(v) is exact if G(v) is a polynomial in v 
of order AN — 1 or less. For our purposes here, G(v) in Eq. (|A1|) is either the -F s (i>, Q 2 ) or the G(v, Q 2 ) on the l.h.s. 
of Eq. ([50)1 or Eq. (|3"Tj) . whereas <?(s) in Eq. (|A1|) is the surrogate for either fc//(s, r)/ (s) + kf g (s,T)g {s) found in 
the r.h.s. of Eq. ([30]) or kgg(s,T)go(s) + (s, t)/o(s) found in the r.h.s. of Eq. (|3ip . Since we must evaluate <?(s) &t 
complex values of s, this necessarily implies that we must evaluate /o(s) and <?o(s) — the Laplace transforms of F s q(v) 
and Gq(v), respectively — at complex values of s. As shown in Ref. to insure numerical accuracy we must be 
able to evaluate g(s) in Eq. (|A1|) to arbitrary accuracy. Thus we must know the Laplace transforms /o(s) and go{s) 
analytically and not just as numerical integrations of the form F s o(v)e~ vs dv. The fc's, the coefficient functions 
needed, are known analytically; the potential problem is with fo(s) an d go(s), the starting functions in Laplace space 
s. 

The starting distributions functions normally used are not of the type that have analytic Laplace transforms. To 
get a sufficiently accurate numerical approximation to functions that do have analytic Laplace transforms is again a 
delicate numerical exercise. We found that we could do it sufficiently accurately by using an interpolating polynomial 
of order n=49. Its 50 coefficients were determined by evaluating the original function at 50 points, distributed as the 
zeroes of a 50 th order Chebyshev polynomial, found in the interval (—1, +1) and then linearly transformed to v space 
to lie in the interval 0<u<14(l>a;> 0.83 x 10~ 6 ). These points were chosen to try to minimize the maximum 
interpolation error. We note that even when using Mathematica, caution was needed in order to obtain sufficient 
numerical accuracy with a such a high order polynomial; it had to be evaluated using Horner's method (see Section 
10.14 of Ref. [l6j). since straight forward evaluation of such a high order polynomial will yield numerical nonsense. 

Using 2N = 38 in Eq. ([AT]) , we would have an exact result if either F s (v, Q 2 ) in Eq. (J3(| or G(v, Q 2 ) Eq. ((31]) were 
a polynomial in v of degree 75 or less; see Ref. Q for details. In actual practice, by comparing the results for the 
value of 2N = 38 — the value we used for our numerical evaluations — with very much larger values of 2N that we 
used for estimates of the exact solutions, we found that the fractional accuracy of inversion for both F s (v,Q 2 ) and 
G(v,Q 2 ) was «±lx 10~ n for v > 0.3. Thus, numerical inversion of the Laplace transform in either Eq. ([30]) or 
Eq. (|31j) contributes essentially nothing to our overall error of about ±2 x 10~ 9 , since it's some 2 orders of magnitude 
smaller. We comment that the overall error is essentially completely due to our numerical approximation of the 
starting functions and not the subsequent Laplace transforms of them. Therefore, we could readily reduce this error 
by using more than 50 points in our numerical approximations of the starting distributions, but this would be at the 
expense of more computational time and was felt to be unnecessary. 

A typical time for computing the full x distribution of either F s (x) or G(x) at an arbitrary Q 2 — given the starting 
functions F s q(x) and Gq(x) at Qq — was about 15 seconds, basically proportional to the number of points in x used in 
the numerical approximations of the starting functions and to the number 2N used in the Laplace inversion routine. 
Thus, for most applications, we could easily reduce this time to several seconds, at the expense of some (perhaps 
unneeded) accuracy. The computations in this paper were made on a home PC, a Dell Model Studio XPS435MT, 
using an Intel 2.67 GHz 4 core i7 CPU, running 64 bit Windows Vista, and using Mathematica? [15J in parallel mode. 

For a very fast Mathematica? (.nb) program that accurately calculates all LO MSTW2008 parton distribution 
functions, as well as F^ix) and F s (x) for any Q 2 — using the LO MSTW starting values Q for F s (x), G(x) at Q 2 , = 1 
GeV 2 — send an email request to mblock@northwestern.edu for MSTW. zip. 
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